34th  Annual  Precise  Time  and  Time  Interval  (PTTI)  Meeting 


HILBERT  SPECTRUM  FOR  TIME-DOMAIN 
MEASUREMENT  DATA  AND  ITS  APPLICATION 


H.  M.  Peng  and  P.  C.  Chang 
National  Standard  Time  &  Frequency  Laboratory,  TL, 

12,  Lane  551,  Min-Tsu  Road,  Sec.  5,  Yang-Mei,  Taoyuan,  Taiwan  326 
Tel:  886  3  424  5174;  Fax:  886  3  424  5474;  E-mail:  hmpeng@cht.com.tw 

F.  R.  Chang 

Department  of  Electrical  Engineering,  NTU,  Taiwan 


Abstract 

In  this  paper,  a  new  method  for  analyzing  the  time  domain  data  is  introduced.  As  one 
knows,  the  time  domain  phase  measurements  are  nonstationary  and  the  differencing  technique 
is  usually  adopted  for  generating  stationary  data.  Thus,  the  analysis  tools  for  stationary  data 
are  used  in  such  case.  On  the  other  hand,  some  nonstationary  data  processing  methods  could 
be  used  directly  to  analyze  measurements,  and  time-varying  phenomena  of  data  will  be  revealed 
through  these  methods.  A  new  algorithm  based  on  the  concept  of  the  instantaneous  frequency 
is  described  in  this  paper.  The  strategy  is  described  as  the  following  two  parts.  First,  the 
measurements  are  decomposed  into  components  with  meaningful  instantaneous  frequencies. 
For  this  purpose,  some  restriction  conditions  for  decomposition  from  the  local  viewpoint  are 
necessary.  After  the  process,  the  Hilbert  transform  will  be  applied  to  every  component  and  the 
data  will  be  expressed  in  phasor  form.  Under  this  circumstance,  we  can  represent  the 
amplitude  and  the  instantaneous  frequency  as  functions  of  time  in  a  three-dimensional  plot. 
This  frequency-time  distribution  of  the  amplitude  is  designated  as  the  Hilbert  amplitude 
spectrum,  and  the  squared  values  of  amplitude  could  be  substituted  for  energy  to  produce  the 
Hilbert  energy  spectrum  just  as  well.  To  verify  our  algorithm,  some  simulations  are  performed. 
Based  on  the  simulation  results,  the  time-varying  spectrum  characteristic  can  be  easily 
identified.  Furthermore,  different  types  of  noise  also  can  be  discriminated  by  its  Hilbert 
spectrum.  Besides,  the  marginal  spectrum  of  the  Hilbert  spectrum  offers  a  measure  of  total 
amplitude  or  energy  contribution  from  each  frequency  value,  which  is  more  meaningful  than 
the  Fourier  spectrum  for  nonstationary  data.  Based  on  the  spectrum  information,  the  optimal 
weighting  factors  between  different  measurements  will  be  further  investigated. 


INTRODUCTION 

For  extracting  the  meaningful  information  and  statistics  from  the  experimental  data,  an  appropriate  data 
analysis  is  necessary.  Generally  speaking,  two  main  purposes  exist  for  data  analysis.  The  first  one  is  to 
determine  the  parameters  needed  to  construct  the  necessary  model,  and  the  other  is  to  confirm  the  model 
constructed  to  represent  the  phenomenon.  Unfortunately,  the  experimental  data  are  always  imperfect,  e.g. 
the  total  data  span  is  too  short,  the  data  are  nonstationary,  or  they  represent  nonlinear  processes.  Facing 
such  data,  we  have  limited  options  to  use  in  the  analysis.  For  example,  the  time-domain  phase 
measurements  are  nonstationary  and  the  differencing  technique  is  usually  adopted  for  generating 
stationary  data.  Thus,  the  analysis  tools  for  stationary  data  are  used  in  such  a  case  [1,2].  On  the  other 
hand,  some  nonstationary  data  processing  methods  could  be  used  directly  to  analyze  measurements  and 
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time-varying  phenomena  of  data  will  be  revealed  through  these  methods. 

In  this  paper,  we  will  use  the  data  analysis  method  introduced  in  [3].  This  scheme  can  be  divided  into 
two  parts.  In  the  first  step,  the  experimental  data  are  decomposed  into  a  collection  of  intrinsic  mode 
functions  (IMF).  This  decomposition  is  viewed  as  an  expansion  of  the  data  in  terms  of  the  IMFs.  In 
other  word,  these  IMFs  are  regarded  as  the  basis  of  that  expansion  which  can  be  linear  or  nonlinear  as 
dictated  by  the  data.  Since  the  IMFs  have  well-behaved  Hilbert  transforms,  the  corresponding 
instantaneous  frequencies  are  calculated.  Thus,  in  the  next  step,  we  could  localize  any  event  on  the  time 
as  well  as  the  frequency  axis.  The  local  energy  and  the  instantaneous  frequency  derived  from  the  IMFs 
give  us  a  full  energy-ffequency-time  distribution  of  the  data,  and  such  a  representation  is  designated  as 
the  Hilbert  spectrum. 

To  verify  the  scheme,  some  simulations  were  conducted.  In  the  first  simulation,  the  white  phase 
modulation  and  white  frequency  modulation  noises  were  tested.  The  marginal  power  spectrum  was 
generated  and  compared  with  the  Fourier  spectrum.  In  the  second  analysis,  these  two  noise  data  were 
concatenated  to  generate  the  time-varying  noise.  Under  this  circumstance,  the  algorithm  was  still 
effective  and  the  spectrum  was  able  to  show  the  time-varying  characteristic. 

This  paper  is  organized  as  follows.  In  the  following  section,  the  intrinsic  mode  function  (IMF)  and  the 
empirical  mode  decomposition  are  described.  Two  types  of  signals  are  used  to  demonstrate  the  results  of 
decomposition.  The  Hilbert  spectrum  and  marginal  average  power  spectrum  are  introduced,  too.  Two 
simulations  are  conducted  and  the  power  spectrum  is  generated  in  the  next  section  to  verify  the  algorithm. 
The  conclusion  is  given  in  the  final  section. 


SUMMARY  OF  THE  EMPIRICAL  MODE  DECOMPOSITION 
METHOD 

Before  introducing  the  empirical  mode  decomposition  proposed  in  [3],  the  concept  of  instantaneous 
frequency  must  be  described  first.  For  an  arbitrary  time  series,  x(t) ,  we  have  its  Hilbert  Transform  [4], 
y(t) ,  as 

y(t)  =  -P[^dt',  (1) 

n  Jo°t-t 

where  P  indicates  the  Cauchy  principal  value.  This  transform  exists  for  all  functions  of  class  Lp .  With 
this  definition,  x(t)  and  y(t)  can  be  used  to  define  an  analytic  signal  z(t) ,  as 


z(t)  =  x{t)  +  iy(t )  =  a(t)e'm , 

in  which 

a(t)  =  [x2  ( t )  +  v 2  {t)Y2 ,  0(t)  =  arctan(^-^)  . 

x(t) 


(2) 

(3) 


In  equation  (2),  the  polar  coordinate  expression  further  clarifies  the  local  nature  of  this  representation:  it 
is  the  best  local  fit  of  an  amplitude  and  phase-varying  trigonometric  function  to  x(t)  .  Furthermore,  the 
instantaneous  frequency  of  x(t)  is  defined  as 


(0  = 


dO(t ) 
dt 


(4) 


Since  the  definition  of  the  instantaneous  frequency  given  in  (4)  is  a  single-valued  function  of  time,  some 
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limitations  on  the  data  are  necessary.  Based  the  discussion  in  [3],  the  authors  proposed  a  class  of 
functions  designated  as  intrinsic  mode  functions  for  the  instantaneous  frequency  to  make  sense.  An 
intrinsic  mode  function  is  a  function  that  satisfies  two  conditions.  The  first  one  is  the  number  of  extrema 
and  the  number  of  zero  crossings  must  either  equal  or  differ  by  one  in  the  whole  data  set.  The  other  one 
is  that  the  mean  value  of  the  envelope  defined  by  the  local  maxima  and  the  envelope  defined  by  local 
minima  is  zero  at  any  point.  Based  on  this  definition  of  IMF,  not  only  are  the  traditional  narrow  band 
requirements  satisfactory,  but  the  unwanted  fluctuations  induced  by  asymmetric  wave  forms  are  reduced. 

For  a  complicated  data  set,  the  decomposition  into  IMF  components  with  meaningful  instantaneous 
frequencies  is  necessary.  This  decomposition  is  based  on  three  assumptions.  First,  the  signal  has  at  least 
two  extrema  -  one  maximum  and  one  minimum.  Next,  the  characteristic  time  scale  is  defined  by  the  time 
lapse  between  the  extrema.  Finally,  if  the  data  are  totally  devoid  of  extrema,  but  contain  only  inflection 
points,  then  it  can  be  differentiated  once  or  more  times  to  reveal  the  extrema.  Based  on  the  discussion  in 
[3],  the  decomposition  can  be  concluded  with  the  following  steps: 

a.  Identify  the  extrema  of  the  data  set  x(t) ,  and  form  the  envelopes  defined  by  the  local 
maxima  and  minima  respectively  by  the  cubic  spline  method. 

b.  Form  the  mean  values  mx(t)  by  averaging  the  upper  envelope  and  lower  envelope,  and 
make  the  differences  between  the  data  and  the  mean  values  to  get  the  first  component 
hx(t)  =  x(t)  - 1 nx  ( t ) . 

c.  If  the  first  component  is  not  an  IMF,  let  hx  ( t )  be  the  new  data  set.  Continue  the  steps  a  and 
b  until  the  first  component  is  an  IMF. 

d.  The  first  IMF  component  is  called  as  cx  ( t ) .  Let  r,  (t)  =  x(t)  -  c,  ( t ) .  Continue  the  steps  a-c 
until  rn  ( t )  is  smaller  than  a  predetermined  value  or  becomes  a  monotonic  function  that  no 
more  IMF  can  be  extracted. 


Based  the  above  algorithm,  two  signals  are  used  to  demonstrate  the  decomposition  results.  The  first 
signal  is  constructed  by  adding  two  cosine  waves  as  5;  ( t )  =  0.2cos(10f)  +  cos(20f) .  This  signal  is  shown 
in  the  left-hand  side  of  Figure  1 .  The  corresponding  IMF  components  cx ,  c2,  and  the  residue  r2  are 
shown  in  the  right-hand  side  of  Figure  1.  Figure  2  shows  the  corresponding  instantaneous  frequencies 
and  amplitudes  of  these  IMFs.  The  second  signal  is  a  time-varying  signal  and  it  is  designated 
Jcos(lOt)  0<t<2 
|cos(20t)  2  <  t  <  4 


as  s2(t)  =  ■ 


Figures  3  and  4  show  the  signal,  IMFs,  and  the  corresponding 


instantaneous  frequencies  and  amplitudes.  From  the  results,  it  is  obvious  that  the  characteristics  of  the 
signals  can  be  extracted  through  the  empirical  mode  decomposition. 


THE  HILBERT  SPECTRUM 

After  the  decomposition  step,  the  data  are  reduced  to  several  IMF  components.  By  performing  the  Hilbert 
transform  on  each  IMF  components,  the  data  can  be  expressed  as  follows: 

x(t)  =  '^aj(t) exp(/ j  (Oj  (t)dt) .  (5) 

7=1 

Here  the  residue,  r  ,  is  omitted  because  it  is  either  a  monotonic  function  or  it  might  be  smaller  than  the 
predetermined  threshold.  Comparing  to  the  Fourier  representation,  it  is  obvious  the  equation  (5)  is  a 
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generalized  Fourier  expansion.  The  time-varying  characteristic  enables  us  to  accommodate  non¬ 
stationary  data.  Furthermore,  equation  (5)  enables  us  to  represent  the  amplitude  and  the  instantaneous 
frequency  as  functions  of  time.  In  other  words,  the  amplitude  can  be  contoured  on  the  frequency-time 
plane.  This  frequency-time  distribution  of  the  amplitude  is  designated  as  the  Hilbert  amplitude  spectrum, 
H(a>,  t ) ,  or  simply  Hilbert  spectrum. 

We  also  can  define  the  marginal  average  power  spectrum  h{(d)  as 


h(co)  =  j^H2(co,t)dt.  (6) 

This  marginal  spectrum  can  be  used  as  a  measure  of  average  power  contribution  from  each  frequency 
value.  It  is  noted  that  the  frequency  in  the  Hilbert  spectrum  has  a  totally  different  meaning  from  Fourier 
spectrum  analysis.  In  the  Fourier  representation,  a  frequency  component  means  a  sine  or  cosine  wave 
persisted  in  the  whole  data  set.  While  the  data  are  nonstationary,  the  Fourier  spectrum  is  meaningless 
physically.  On  the  other  hand,  the  frequency  in  the  Hilbert  spectrum  or  marginal  spectrum  indicates  that 
an  oscillation  with  such  a  frequency  exists.  The  representation  is  still  valid  while  the  data  set  is 
nonstationary. 


THE  SIMULATIONS 

To  verify  the  method,  two  simulations  were  conducted.  In  the  first  simulation,  two  types  of  noises  were 
generated.  The  first  one  was  the  white  phase  modulation  noise  and  the  other  one  was  white  frequency 
modulation  noise  with  the  same  noise  parameter  <7r  =1.0x10“"  .  The  sampling  time  is  1  second  and  the 

time  span  is  400  seconds.  Figures  5  and  7  show  the  noise  plots  and  the  IMF  components,  respectively.  In 
this  case,  the  white  phase  modulation  noise  can  be  decomposed  into  nine  IMF  components  and  the  white 
frequency  phase  modulation  noise  has  seven  IMFs.  Using  equation  (6),  we  can  generate  the  Hilbert 
marginal  power  spectrum  and  the  results  are  shown  in  Figures  6  and  8,  respectively.  The  corresponding 
Fourier  power  spectra  are  also  shown  in  the  same  figures  for  comparison. 

In  the  second  simulation,  we  concatenate  the  two  noises  as  a  new  time-varying  noise.  The  IMF 
components  are  shown  in  Figure  9  and  the  Hilbert  power  spectrum  is  shown  in  Figure  10.  From  Figure 
10,  we  can  observe  that  the  power  density  is  increasing  at  low  frequency  with  time.  This  feature  can  be 
used  to  demonstrate  that  the  time-varying  characteristics  of  the  concatenated  noise  are  revealed  via  the 
decomposition  method. 


CONCLUSION 

In  this  paper,  we  use  the  analysis  method  proposed  in  [3]  to  analyze  the  time  domain  measurements.  The 
data  can  be  decomposed  into  several  well-behaved  functions  that  have  meaningful  instantaneous 
frequencies.  Based  on  the  decomposition,  we  can  observe  the  frequency  characteristics  of  the 
experimental  data.  Moreover,  even  if  the  frequency  components  of  the  data  are  time-varying,  the  analysis 
scheme  is  still  valid.  By  using  this  tool,  we  can  observe  the  local  phenomena  of  the  measurements. 

REFERENCES 

[1]  D.  B.  Sullivan,  D.  W.  Allan,  D.  A.  Howe,  and  F.  L.  Walls  (editors),  1990,  Characterization  of 
Clocks  and  Oscillators,  NIST  Technical  Note  1337  (National  Institute  of  Standards  and 


504 


34,h  Annual  Precise  Time  and  Time  Interval  (PTTI)  Meeting 


Technology). 

[2]  J.  Rutman  and  F.  L.  Walls,  1991,  “Characterization  of  frequency  stability >  in  precision  frequency 
sources,  ”  Proceedings  of  the  IEEE,  79,  952-960. 

[3]  N.  E.  Huang  et  al.,  1998,  “The  empirical  mode  decomposition  and  the  Hilbert  spectrum  for 
nonlinear  and  non- stationary  time  series  analysis,  ”  Proceedings  of  the  Royal  Society  of  London, 
454,  903-995. 

[4]  J.  S.  Bendat  and  A.  G.  Piersol,  1986,  Random  Data:  Analysis  and  Measurement  Procedures, 
second  edition  (John  Wiley  &  Sons,  New  York). 


s,(t)  0 


3  4 


decompositoin 


Figure  1.  s,(t)  and  its  corresponding  IMF  components. 
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Figure  2.  The  instantaneous  frequencies  and  amplitudes  of  IMFs  in  sx(t) . 
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Figure  3.  s2(t)  and  its  corresponding  IMF  components. 
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Figure  4.  The  instantaneous  frequencies  and  amplitudes  of  IMFs  in  s2(t) . 
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Figure  5.  White  phase  modulation  noise  and  its  IMF  components. 
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Fourier  power  psectrum  (red)  and  Hilbert  marginal  power  spectrum  (blue) 
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Figure  6.  Marginal  spectrum  and  power  spectrum  for  white  phase  modulation  noise. 
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Figure  7.  White  frequency  modulation  noise  and  its  IMF  components. 
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Fourier  power  spectrum  (red)  and  Hilbert  marginal  power  spectrum  (blue) 
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Figure  8.  Marginal  spectrum  and  power  spectrum  for  white  frequency  modulation  noise. 
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Figure  9.  IMF  components  of  the  concatenated  noise. 
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Figure  10.  Flilbert  power  spectrum  for  the  concatenated  noise. 
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